#*************        丁香园-NHANES     *************
#*************        Analyze Code      *************
#*************       分析结果复现       *************
#*************                          *************

# Paper01-美国成年人膳食中类胡萝素与认知功能的关系, NHANES 2011-2014
# Paper01-Dietary carotenoids and cognitive function among US adults, NHANES 2011–2014

#### 0.准备好环境 ####
library(gtsummary)
library(survey)
library(haven)
library(tableone)
library(plyr)
library(dplyr) # 链接：https://dplyr.tidyverse.org/reference/mutate.html
library(tidyverse)
library(arsenal) 
setwd("G:/BaiduNetdiskDownload/NHANES")
#setwd("J:/nhanes/数据分析/NHANES_20221011") #需要转换为自己的数据读取路径

#### 一. 定位数据模块和变量，获取源数据 ####
##### 1.1 DEMO-人口学数据提取 #####
### 1.1 提取 Component文件
# NHANES官网链接：https://wwwn.cdc.gov/Nchs/Nhanes/2011-2012/DEMO_G.htm
demo.g <- read_xpt("2011-2012/Demographics/demo_g.xpt")#参见上述设置默认路径，文件名称后缀不同"g,h"
demo.h <- read_xpt("2013-2014/Demographics/demo_h.xpt")#参见上述设置默认路径，文件名称后缀不同"g,h"

### 1.2 提取研究所需要的变量
# 年龄-RIDAGEYR; 性别-RIAGENDR; 种族-RIDRETH3; 教育程度-DMDEDUC2; 贫困程度-INDFMPIR;
demo.data.file <- dplyr::bind_rows(list(demo.g, demo.h))
demo.data <- demo.data.file[,c('SEQN', 'RIDAGEYR', 'RIAGENDR', 'RIDRETH3', 'DMDEDUC2', 'INDFMPIR')]

##### 1.2 SMQ-吸烟数据提取 #####
### 1.1 提取 Component文件
smq.g <- read_xpt("2011-2012/Questionnaire/smq_g.xpt")
smq.h <- read_xpt("2013-2014/Questionnaire/smq_h.xpt")

### 2. 提取研究所需要的变量
# 是否吸烟至少100支-SMQ020; 现在是否吸烟-SMQ040; 多久前开始戒烟-SMQ050Q;
# 多久前开始戒烟的时间单位（天、周、月、年）-SMQ050U;
smq.data.file <- dplyr::bind_rows(list(smq.g, smq.h))
smq.data <- smq.data.file[,c('SEQN', 'SMQ020', 'SMQ040', 'SMQ050Q', 'SMQ050U')]

##### 1.3 ALQ-饮酒数据提取 ######
### 1.提取 Component文件
alq.g <- read_xpt("2011-2012/Questionnaire/alq_g.xpt")
alq.h <- read_xpt("2013-2014/Questionnaire/alq_h.xpt")

### 2.提取研究所需要的变量
# 每年至少喝12杯酒-ALQ101
# 一生中至少喝过12次酒-ALQ110；
# 过去12个月多久喝一次酒-ALQ120Q
# 过去12个月饮酒频率的单位（周，月，年）-ALQ120U;

alq.data.file <- dplyr::bind_rows(list(alq.g, alq.h))
alq.data <- alq.data.file[,c('SEQN', 'ALQ101', 'ALQ110', 'ALQ120Q', 'ALQ120U')]

##### 1.4 BMX-BMI、腰围数据提取 #####
### 1.提取 Component文件
bmx.g <- read_xpt("2011-2012/Examination/bmx_g.xpt") #注意是 Examination 的类别
bmx.h <- read_xpt("2013-2014/Examination/bmx_h.xpt") #注意是 Examination 的类别

### 2.提取研究所需要的变量 BMI-BMXBMI; 腰围-BMXWAIST
bmx.data.file <- dplyr::bind_rows(list(bmx.g, bmx.h))
bmx.data <- bmx.data.file[,c('SEQN', 'BMXBMI', 'BMXWAIST')]

##### 1.5 饮食-L&Z 以及总卡路里数据提取(第1天 & 第2天) #####
# 第1天
dr1tot.g <- read_xpt('2011-2012/Dietary/dr1tot_g.xpt')
dr1tot.h <- read_xpt('2013-2014/Dietary/dr1tot_h.xpt')

dr1tot.data.file <- dplyr::bind_rows(list(dr1tot.g, dr1tot.h))
dr1tot.data <- dr1tot.data.file[,c('SEQN', 'DR1TLZ', 'DR1TKCAL')]

# 第2天
dr2tot.g <- read_xpt('2011-2012/Dietary/dr2tot_g.xpt')
dr2tot.h <- read_xpt('2013-2014/Dietary/dr2tot_h.xpt')

dr2tot.data.file <- dplyr::bind_rows(list(dr2tot.g, dr2tot.h))
dr2tot.data <- dr2tot.data.file[,c('SEQN', 'DR2TLZ', 'DR2TKCAL')]

# 合并第1天&第2天的变量
dr.data <- merge(dr2tot.data, dr1tot.data)

##### 1.6 膳食补充剂-L&Z数据提取(第1天 & 第2天) ######
# 第1天
ds1tot.g <- read_xpt('2011-2012/Dietary/ds1tot_g.xpt')
ds1tot.h <- read_xpt('2013-2014/Dietary/ds1tot_h.xpt')

ds1tot.data.file <- dplyr::bind_rows(list(ds1tot.g, ds1tot.h))
ds1tot.data <- ds1tot.data.file[,c('SEQN', 'DS1TLZ')]

# 第2天
ds2tot.g <- read_xpt('2011-2012/Dietary/ds2tot_g.xpt')
ds2tot.h <- read_xpt('2013-2014/Dietary/ds2tot_h.xpt')

ds2tot.data.file <- dplyr::bind_rows(list(ds2tot.g, ds2tot.h))
ds2tot.data <- ds2tot.data.file[,c('SEQN', 'DS2TLZ')]

# 合并第1天&第2天的变量
ds.data <- merge(ds1tot.data, ds2tot.data)
# View(ds.data)

##### 1.7 CFQ-认知数据提取 #####
### 1.提取 Component文件
cfq.g <- read_xpt('2011-2012/Questionnaire/cfq_g.xpt')
cfq.h <- read_xpt('2013-2014/Questionnaire/cfq_h.xpt')

### 2.提取研究所需要的变量
#CERAD 词汇第1次即刻记忆测验	CFDCST1
#CERAD 词汇第2次即刻记忆测验	CFDCST2
#CERAD 词汇第3次即刻记忆测验	CFDCST3
#CERAD 3次词汇即刻记忆测验总分	CERAD_1+CERAD_2+CERAD_3
#CERAD 延迟回忆评分	CFDCSR
#语言流畅性评分Animal Fluency	CFDAST
#数字符号替代测试(DSST)	CFDDS

cfq.data.file <- dplyr::bind_rows(list(cfq.g, cfq.h))
cfq.data <- cfq.data.file[,c('SEQN', 'CFDCST1', 'CFDCST2', 'CFDCST3',
                             'CFDCSR', 'CFDAST', 'CFDDS')]

#### 2 提取分析相关变量（权重等，暂时为复现paper结果而提取） ####
##### 2.1 权重变量 ##### 
# 找到权重变量
# DEMO->MEC-Dietary->MEC-CFQ
# 饮食中的权重：https://wwwn.cdc.gov/nchs/nhanes/2013-2014/DR1TOT_H.htm

# SMQ: https://wwwn.cdc.gov/Nchs/Nhanes/2013-2014/SMQ_H.htm
# SMQ 中没有权重信息：The NHANES full sample 2-Year MEC Exam Weights (WTMEC2YR) should be used to analyze the 2013-14 SMQ variables in conjunction with the laboratory measurements on tobacco exposure or other examination measurements.

# CFQ: https://wwwn.cdc.gov/Nchs/Nhanes/2011-2012/CFQ_G.htm
# CFQ 中没有权重信息：MEC exam sample weights should be used for the analyses 
# and are available on the NHANES website 

weight.data <- dr1tot.data.file[,c('SEQN', 'WTDRD1')]
weight.data$WTDRD1 <- weight.data$WTDRD1/2 #权重的计算方式，需要按 Cycle 进行平均

##### 2.2 复杂抽样的其他变量-DEMO #####
survey.design.data <- demo.data.file[,c('SEQN', 'SDMVPSU', 'SDMVSTRA')]

#### 3.合并上述所有数据（把列拼接起来）-Output #### 
##### 3.1 合并步骤1 & 2中提取的数据 #####
output <- plyr::join_all(list(demo.data, smq.data, alq.data, bmx.data,
                              dr.data, ds.data, cfq.data, weight.data, survey.design.data),
                         by='SEQN', type='full')

##### 3.2 针对合并后的数据进行质控，确定人数能对得上 #####
# 1.根据 Table1 进行质控
data.age.60 <- subset.data.frame(output, RIDAGEYR >= 60 )
dim(data.age.60) #[1] 3632    ✅

# 2. 确认下权重
weight.non.na.index <- which(!is.na(data.age.60$WTDRD1))
sum(data.age.60$WTDRD1[weight.non.na.index]) # 计算权重-59659564 ✅


#### 4.根据文章的筛选策略进行初步筛选（仅根据原始数据） ####
# 原始筛选策略：The analyses presented here include a pooled sample of all individuals aged
# 60 years or older in each of the two survey cycles, with completed cognitive performance
# test results, dietary L and Z intake information, and information on important
# con-founders(age, sex, race/ethnicity, smoking, educational attainment).

# 注：这里不需要引号，e.g. 直接写 CFDCSR 就可以了（视频中存在错误）-12.06更新
paper.data <- subset.data.frame(output, RIDAGEYR >= 60 &
                                  (!is.na(CFDCSR)) & #CERAD 延迟回忆评分
                                  (!is.na(CFDAST)) & #语言流畅性评分Animal Fluency
                                  (!is.na(CFDDS)) & #数字符号替代测试(DSST)
                                  (!is.na(DR1TLZ))& #饮食-叶黄素 & 玉米黄质摄取
                                  (!is.na(DR1TKCAL))& #能量摄入总量-Energy (kcal)
                                  (!is.na(RIDAGEYR)) & #年龄
                                  (!is.na(RIAGENDR)) & #性别
                                  (!is.na(RIDRETH3)) & #种族
                                  (!is.na(DMDEDUC2)) & #教育程度
                                  (!is.na(SMQ020)) #是否吸烟至少100支
)
dim(paper.data)

#### 5.衍生得到相应的变量，同步衍生 ####
##### 5.1 性别-Sex#####
paper.data$Sex <- ifelse(paper.data$RIAGENDR == 1, 'male', 'female')
##### 5.2 年龄-age.group #####
paper.data$age.group <- ifelse(paper.data$RIDAGEYR >= 60 & paper.data$RIDAGEYR < 69, '60-69 years',
                               ifelse(paper.data$RIDAGEYR >=70 & paper.data$RIDAGEYR < 79, '70-79 years',
                                      '80+ years'))
##### 5.3 种族-race #####
race <- recode_factor(paper.data$RIDRETH3, 
                      `1` = 'Mexican American',
                      `2` = 'Other Hispanic',
                      `3` = 'Non-Hispanic White',
                      `4` = 'Non-Hispanic Black',
                      `6` = 'Other/multiracial',
                      `7` = 'Other/multiracial'
)
paper.data$race <- race

##### 5.4 复杂转换-饮酒 alq.group #####
# 仅衍生有无饮酒
# 区分男性和女性不同的饮酒标准，drinks per day (dkspd) ：
#### 是否酗酒：https://bmcgastroenterol.biomedcentral.com/articles/10.1186/s12876-022-02430-7
#### 涉及男性和女性的饮酒标准不一致：https://www.nature.com/articles/s41598-017-02426-4  


# 衍生方式：Paper01-美国成年人膳食中类胡萝素与认知功能的关系, NHANES 2011-2014
# 衍生结果为：Non-drinker, 1-5 drinks/month, 5-10 drinks/month, 10+ drinks/month
# https://wwwn.cdc.gov/Nchs/Nhanes/2011-2012/ALQ_G.htm#ALQ120U （课程使用的paper）
# ALQ120Q, 具体的数值;  777:refused, 999:don't know
# ALQ120U, 1:week, 2:month, 3:year, 7:refused, 9:don't know


# 统一单位为月
# week-month: *4; year-month:/12
ori.alq.unit <- paper.data$ALQ120U

trans.unit.month <- ifelse(ori.alq.unit == 1, 4, 
                           ifelse(ori.alq.unit == 3, 1/12, 
                                  ifelse(ori.alq.unit == 7|ori.alq.unit == 9, NA, 1)))
paper.data$trans.unit.month <- trans.unit.month
# View(paper.data)

# 统一数值为月饮酒量-quantity
ori.alq.quantity <- paper.data$ALQ120Q
trans.quantity.month <- ifelse(ori.alq.quantity >= 0,  
                               ori.alq.quantity * trans.unit.month, NA)
paper.data$trans.quantity.month <- trans.quantity.month


# 将饮酒量的结果划分为4档：Non-drinker, 1-5 drinks/month, 5-10 drinks/month, 10+ drinks/month
# 首先根据转化得到的trans.quantity.month，将饮酒量划分3档
# 利用 paper.data$ALQ101
alq101 <- paper.data$ALQ101
trans.quantity.month.factor <- ifelse(trans.quantity.month >=1 & trans.quantity.month <5, '1-5 drinks/month',
                                      ifelse(trans.quantity.month >=5 & trans.quantity.month <10, '5-10 drinks/month',
                                             ifelse(trans.quantity.month >=10, '10+ drinks/month', 'wait')))

paper.data$trans.quantity.month.factor <- trans.quantity.month.factor

ddply(paper.data, .(ALQ101, trans.quantity.month.factor), summarise, n = length(SEQN))

# ALQ101 回答是1，但是 trans.quantity.month.factor 没有具体的值，则将 trans.quantity.month.factor 取值为1
index.1 <- which((trans.quantity.month.factor=='wait' | is.na(trans.quantity.month.factor)) & alq101 == 1) #844
trans.quantity.month.factor[index.1] <- '1-5 drinks/month'

paper.data$trans.quantity.month.factor <- trans.quantity.month.factor

index.less.1 <- which(alq101 == 2)
trans.quantity.month.factor[index.less.1] <- 'Non-drinker'

table(trans.quantity.month.factor)
# 根据是否饮酒补充 Non-drinker
paper.data$alq.group <- trans.quantity.month.factor

##### 5.5 Total 卡路里-total.calories #####
# day1、day2 都有，则取2天平均，否则，取第1天的值
# View(paper.data[,c('DR1TKCAL', 'DR2TKCAL')])
# 先计算2天的平均值，其中第二天为 NA 的，平均值也是 NA
total.calories <- apply(paper.data[,c('DR1TKCAL', 'DR2TKCAL')], 1, mean)
paper.data$total.calories <- total.calories

# 把第二天为 NA 的值计算出来，用第一天的值作为平均值
day.2.na.index <- which(is.na(paper.data$DR2TKCAL))
total.calories[day.2.na.index] <- paper.data$DR1TKCAL[day.2.na.index]

# 添加到paper.data中，再次确认
paper.data$total.calories <- total.calories


##### 5.6 Total 饮食中 L&Z-total.dr.LZ#####
# https://wwwn.cdc.gov/Nchs/Nhanes/2011-2012/DR2TOT_G.htm#DR2TLZ
# day1、day2 都有，则取2天平均，否则，取第一天的值
# 先计算2天的平均值，其中第二天为 NA 的，平均值也是 NA
total.dr.LZ <- apply(paper.data[,c('DR1TLZ', 'DR2TLZ')], 1, mean)
paper.data$total.dr.LZ <- total.dr.LZ

# 把第二天为 NA 的值计算出来，用第一天的值作为平均值
day.2.na.index <- which(is.na(paper.data$total.dr.LZ))
total.dr.LZ[day.2.na.index] <- paper.data$DR1TLZ[day.2.na.index]

# 添加到paper.data中，再次确认
paper.data$total.dr.LZ <- total.dr.LZ

##### 5.7 认知评分 #####
#CERAD 词汇第1次即刻记忆测验	CFDCST1
#CERAD 词汇第2次即刻记忆测验	CFDCST2
#CERAD 词汇第3次即刻记忆测验	CFDCST3
#CERAD 3次词汇即刻记忆测验总分	CERAD_1+CERAD_2+CERAD_3

CERAD.total <- apply(paper.data[,c('CFDCST1', 'CFDCST2', 'CFDCST3')], 1, sum)
paper.data$CERAD.total <- CERAD.total

##### 5.8 BMI.group #####
paper.data$BMI.group <- ifelse(paper.data$BMXBMI <18.5, 'Underweight(<18.5)',
                               ifelse(paper.data$BMXBMI >=18.5 & paper.data$BMXBMI < 25, 'Normal(18.5 to <25)',
                                      ifelse(paper.data$BMXBMI >=25 & paper.data$BMXBMI < 30, 'Overweight(25 to <30)',
                                             'Obese(30 or greater)')))

##### 5.9 smoke.group #####
# 本次衍生采取常规简单的定义
# SMQ020:Smoked at least 100 cigarettes in life, 1:Yes, 2:No
# SMQ040:Do you now smoke cigarettes, 1:Every day, 2:Some days	

# 使用 dplyr 的 mutate 函数可以 Create, modify, and delete 列，下面使用 mutate 衍生 smoke group 变量
# 波浪号 ~ 前面是条件，后面是满足这个条件的取值
paper.data <- mutate(paper.data, smoke.group = case_when(
  SMQ020 == 2 ~ 'Never smoker',
  SMQ020 == 1 & SMQ040 == 3 ~ 'Former smoker',
  SMQ020 == 1 & SMQ040 <= 2 ~ 'Current smoker'
))


##### 5.10 education.attainment #####
# 仅仅是单个变量的转换，可以使用 recode_factor，注意，原先的取值是数值或者中文的结果，要使用``标注
# 要注意 7\9\77\99这类的取值，不处理的话直接是 NA, warning 提示：Unreplaced values treated as NA as `.x` is not compatible.

education.attainment <- recode_factor(paper.data$DMDEDUC2, 
                                      `1` = 'Less Than 9th Grade',
                                      `2` = '9-11th Grade',
                                      `3`= 'High School Grad/GED',
                                      `4`= 'Some College or AA degree',
                                      `5`= 'College Graduate or above')

paper.data$education.attainment <- education.attainment


##### 5.11 衍生 L & Z 中剩余的 2 个连续变量 #####
# 膳食补充剂中摄入 L&Z 的摄入量，supplement L and Z intake，命名为 total.ds.LZ，进一步计算出总体 L & Z 的摄入量，命名为 total.dr.ds.LZ
##### 5.11.1 计算total.ds.LZ #####
total.ds.LZ <- apply(paper.data[,c('DS1TLZ', 'DS2TLZ')], 1, mean)
paper.data$total.ds.LZ <- total.ds.LZ

# 把第二天为 NA 的值计算出来，用第一天的值作为平均值
day.2.na.index <- which(is.na(paper.data$total.ds.LZ))
total.ds.LZ[day.2.na.index] <- paper.data$DS1TLZ[day.2.na.index]

# 添加到paper.data中
paper.data$total.ds.LZ <- total.ds.LZ

##### 5.11.2 计算total.dr.ds.LZ #####
# 饮食和膳食补充剂加总
total.dr.ds.LZ <- apply(paper.data[,c('total.dr.LZ', 'total.ds.LZ')], 1, sum)
paper.data$total.dr.ds.LZ <- total.dr.ds.LZ

# Review 一下结果，这一步很重要，在每一步新生成变量的时候，都要看一下原始的结果
# View(paper.data[,c('total.dr.LZ', 'total.ds.LZ', 'total.dr.ds.LZ')]) 

# 将 total为 NA 的结果提取出来，如果一个人仅有 DR，加总的结果还是 NA 因此要另外处理
ds.dr.LZ.na.index <- which(is.na(paper.data$total.dr.ds.LZ))
total.dr.ds.LZ[ds.dr.LZ.na.index] <- paper.data$total.dr.LZ[ds.dr.LZ.na.index]

# 将更新后的total.dr.ds.LZ加入到 paper.data 的结果中
paper.data$total.dr.ds.LZ <- total.dr.ds.LZ

# Review 一下结果，这一步很重要，在每一步新生成变量的时候，都要看一下原始的结果
# View(paper.data[,c('total.dr.LZ', 'total.ds.LZ', 'total.dr.ds.LZ')]) 


# 注意，我们之前还差一步，就是将 LZ 摄入的单位从 mcg 调整为 mg
paper.data$total.dr.LZ <- paper.data$total.dr.LZ/1000
paper.data$total.ds.LZ <- paper.data$total.ds.LZ/1000
paper.data$total.dr.ds.LZ <- paper.data$total.dr.ds.LZ/1000

which(!is.na(paper.data$total.ds.LZ))# 有膳食补充剂数据的一共 471 人


#### 6.再次筛选 ####
# 根据全部的结果，去掉核心变量缺失的行
# there are 2796 participants who had information on cognitive
# performance and dietary L and Z intake (some missing supplement intake), and who had information on important covariates (sex, age, race/ethnicity, smoking, education).
paper.data <- subset.data.frame(paper.data, 
                                (!is.na(total.dr.LZ)) & 
                                  (!is.na(total.dr.ds.LZ)) &
                                  (!is.na(smoke.group))& 
                                  (!is.na(education.attainment))&
                                  (!is.na(CFDCSR)) & #CERAD 延迟回忆评分
                                  (!is.na(CFDAST)) & #语言流畅性评分Animal Fluency
                                  (!is.na(CFDDS))) #数字符号替代测试(DSST)) #是否吸烟至少100支
dim(paper.data) #2712
colnames(paper.data)
# analyze.variable <- c("SEQN", "WTDRD1", "SDMVPSU", "SDMVSTRA",
#                       "Sex", "RIDAGEYR", "age.group", "race", "education.attainment", "INDFMPIR",
#                       "alq.group", "smoke.group","BMXBMI", "BMI.group", "BMXWAIST", 
#                       "total.calories", "total.dr.LZ", "total.ds.LZ", "total.dr.ds.LZ",
#                       "CFDCST1", "CFDCST2", "CFDCST3",  "CERAD.total", 
#                       "CFDCSR", "CFDAST", "CFDDS")
# 
# paper.data <- paper.data[, analyze.variable]
# colnames(paper.data) <- c("SEQN", "WTDRD1", "SDMVPSU", "SDMVSTRA",
#                           "Sex", 'Age', "Age.group", "Race", "Education.attainment", "PIR",
#                           "Alq.group", "Smoke.group", "BMI", "BMI.group", "Waist", 
#                           "Total.calories", "Dietary.LZ", "Supplement.LZ", "Total.LZ",
#                           "CERAD1", "CERAD2", "CERAD3", "CERAD.total",
#                           "CERAD.delay.recall", "Animal.Fluency", "DSST")







#### 二、Paper Table ####
#### 1. 分析数据准备 ####
# 生成复杂抽样的对象
NHANES_design <- svydesign(data = paper.data, ids = ~SDMVPSU, strata = ~SDMVSTRA, nest = TRUE, weights = ~WTDRD1, survey.lonely.psu = "adjust") 

# 修饰表格1-修改分类变量取值出现的顺序

# 按照临床含义排序-1饮酒、吸烟
paper.data$Alq.group <- factor(paper.data$Alq.group, 
                               levels = c('1-5 drinks/month', 
                                          '5-10 drinks/month',
                                          '10+ drinks/month',
                                          'Non-drinker'))

paper.data$BMI.group <- factor(paper.data$BMI.group, 
                               levels = c('Underweight(<18.5)', 
                                          'Normal(18.5 to <25)',
                                          'Overweight(25 to <30)',
                                          'Obese(30 or greater)'))

# 按照取值的顺序-Race
levels(fct_infreq(paper.data$Race)) # 从大到小的顺序
# levels(fct_rev(fct_infreq(paper.data$Race))) # 从小到大的顺序

paper.data$Race <- factor(paper.data$Race, levels(fct_infreq(paper.data$Race)))# 从大到小的顺序
# paper.data$Race <- factor(paper.data$Race, levels(fct_rev(fct_infreq(paper.data$Race))))# 从小到大的顺序

# 加权下的排序
svytotal(~ Race, NHANES_design, na.rm=TRUE) # 离散性变量

NHANES_design <- svydesign(data = paper.data, ids = ~SDMVPSU, strata = ~SDMVSTRA, nest = TRUE, weights = ~WTDRD1, survey.lonely.psu = "adjust") 

#### 2. Table1 ####
tbl_svysummary(NHANES_design,  missing = 'no', 
               include = c(Age.group, Sex,  Race, BMI.group, Alq.group, Smoke.group, Education.attainment,
                           Age, PIR, BMI, Waist, Total.calories, Dietary.LZ, Supplement.LZ, Total.LZ,
                           CERAD1, CERAD2, CERAD3, CERAD.total, CERAD.delay.recall, Animal.Fluency, DSST)) 

tbl_svysummary(NHANES_design,  missing = 'no',
               include = c(Age.group, Sex,  Race, BMI.group, Alq.group, Smoke.group, Education.attainment,
                           Age, PIR, BMI, Waist, Total.calories, Dietary.LZ, Supplement.LZ, Total.LZ,
                           CERAD1, CERAD2, CERAD3, CERAD.total, CERAD.delay.recall, Animal.Fluency, DSST),
               statistic = list(all_continuous()  ~ "{median} ({p25}, {p75})", # 或者常用的："{mean} ({sd})"
                                all_categorical() ~ "{n_unweighted} ({p}%)"))%>% # 替换为非加权的 n
  modify_header(all_stat_cols() ~ "**{level}**, N = {n_unweighted} ({style_percent(p)}%)")


# 也可以用 tableby+weights参数、或者 svymean、svyquantile、svytotal 得到相同的结果
tab1 <- tableby(~ BMI.group + CERAD.total, data = paper.data, weights=WTDRD1,subset = WTDRD1>0) 
summary(tab1, text=TRUE)

svymean(~ CERAD.total, NHANES_design, na.rm=TRUE)
svyquantile(~ CERAD.total, NHANES_design, na.rm=TRUE, quantiles=c(.25,.5,.75), interval.type="mean") # 连续性变量：20, (17, 23)
svytotal(~ BMI.group, NHANES_design, na.rm=TRUE) # 离散性变量

#### 3. Table2 ####
##### 3.1 四分位数变量衍生 ####
# 无权重的情况，对连续变量进行离散分类
paper.data$Dietary.LZ.quantile.var <- cut(paper.data$Dietary.LZ,
                                        breaks = quantile(paper.data$Dietary.LZ),
                                        labels = c('Q1', 'Q2', 'Q3', 'Q4'))
                                           
# View(paper.data[,c('Dietary.LZ.quantile.var', 'Dietary.LZ')])

# 加权的情况，对连续变量进行离散分类
# https://rdrr.io/cran/survey/man/svyquantile.html
Dietary.LZ.quantile.res <- svyquantile(~ Dietary.LZ, NHANES_design, quantiles = c(0, 0.25, 0.5, 0.75, 1))
paper.data$Dietary.LZ.quantile.var <- cut(paper.data$Dietary.LZ,
                                          breaks = Dietary.LZ.quantile.res$Dietary.LZ[,'quantile'],
                                          labels = c('Q1', 'Q2', 'Q3', 'Q4'))
                                           
# View(paper.data[,c('Dietary.LZ.quantile.var', 'Dietary.LZ')])
# forcats::fct_count(paper.data$Dietary.LZ.quantile.var)
paper.data$Dietary.LZ.quantile.var[which(is.na(paper.data$Dietary.LZ.quantile.var))] <- 'Q1'
#  # 将最小值替换为 Q1，以免为 NA

# 加权情况
# 注意：更新来 paper.data 后，需要更新 design 
NHANES_design <- svydesign(data = paper.data, ids = ~SDMVPSU, strata = ~SDMVSTRA, nest = TRUE, weights = ~WTDRD1, survey.lonely.psu = "adjust") 

##### 3.2 四分位数 Table2 一键生成结果 ####
# 使用函数 tbl_svysummary，来一键生成结果
# https://www.danieldsjoberg.com/gtsummary/reference/tbl_svysummary.html

tbl_svysummary(NHANES_design, by = Dietary.LZ.quantile.var, missing = 'no',
               include = c(Dietary.LZ, Animal.Fluency, CERAD1, CERAD2, CERAD3, CERAD.total, CERAD.delay.recall,DSST),
               label = list(Dietary.LZ ~ 'Dietary L and Z intake', Animal.Fluency ~ 'Animal Fluency: Score Total')) %>%     # 设置变量的展示标签          
  modify_header(all_stat_cols() ~ "**{level}**, N = {n_unweighted} ({style_percent(p)}%)") %>%               
  add_p() 

##### 3.3 四分位数 Table2 subset-子群分析 ####
# 以男性为例
NHANES_design %>%
  subset(Sex == 'male') %>%
  tbl_svysummary(by = Dietary.LZ.quantile.var, missing = 'no',
                 include = c(Dietary.LZ),
                 label = list(Dietary.LZ ~ 'Men only')) # 设置变量的展示标签

tbl_svysummary(subset(NHANES_design, Sex == 'male'), 
               by = Dietary.LZ.quantile.var, missing = 'no',
               include = c(Dietary.LZ),
               label = list(Dietary.LZ ~ 'Men only')) # 设置变量的展示标签

                 
# 注意：p value 的计算方式，其检验的方法，一共4种，分别在下述条件下使用：
# Tests default to "kruskal.test" for continuous variables ("wilcox.test" when "by" variable has two levels), 
# "chisq.test.no.correct" for categorical variables with all expected cell counts >=5, 
# and "fisher.test" for categorical variables with any expected cell count <5.

tbl_svysummary(NHANES_design, by = Dietary.LZ.quantile.var, missing = 'no',
               include = c(Dietary.LZ, Animal.Fluency, CERAD1, CERAD2, CERAD3, CERAD.total, CERAD.delay.recall,DSST),
               label = list(Dietary.LZ ~ 'Dietary L and Z intake', Animal.Fluency ~ 'Animal Fluency: Score Total')) %>%     # 设置变量的展示标签          
  modify_header(all_stat_cols() ~ "**{level}**, N = {n_unweighted} ({style_percent(p)}%)") %>%               
  add_p() %>%
  as_flex_table() %>% # 导出 Word
  flextable::save_as_docx(path = 'Table2_Result.docx')

#### 4. Table3  #####
# 使用函数 tbl_regression，来一键生成结果
# https://www.danieldsjoberg.com/gtsummary/articles/tbl_regression.html
# Y1 = Dietary L and Z (mg/day)
# Age-adjusted
m1 <- svyglm(CERAD.delay.recall ~ Dietary.LZ + Age, design = NHANES_design) 
tbl_regression(m1)
tbl_regression(m1, exponentiate = TRUE)
tbl_regression(m1, exponentiate = TRUE, include = c(Dietary.LZ)) # 只呈现 Dietary.LZ

# Fully adjusted
m1 <- svyglm(CERAD.delay.recall ~ Dietary.LZ + Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment, 
             design = NHANES_design)
tbl_regression(m1, include = c(Dietary.LZ)) # 只呈现 Dietary.LZ
tbl_regression(m1, exponentiate = TRUE, include = c(Dietary.LZ)) # 只呈现 Dietary.LZ

# Q4 vs Q1 
m1 <- svyglm(CERAD.delay.recall ~ Dietary.LZ.quantile.var + Age + Sex + BMI + Alq.group + Smoke.group + PIR + Education.attainment, 
             design = NHANES_design)
tbl_regression(m1, exponentiate = TRUE, include = c(Dietary.LZ.quantile.var)) # 只呈现 Dietary.LZ


